# ============================================================================================================
# Replication Do-File for AJPS-Article "Voting against your constituents? How lobbying affects representation"
# by Nathalie Giger and Heike Kl�ver, 03.02.2015
# ============================================================================================================


# Loading Packages
library(foreign)
library(nlme)
library(lattice)
library(Matrix)
library(lme4)
library(car)
library(arm)
library(memisc)
library(mice)
library(Zelig)
library(Amelia)
library(plotrix)      
options(scipen=3)      
library(graphics)

# =============================
# Figure 1
# =============================


# Figure 1a
# =============================
rm(list=ls(all=TRUE))
setwd("C:/Users/ba9ze3-1/Dropbox/Dokumente/Papers/Schweiz IG MP Paper/Final Submission/Giger_Kluever_Replication")
swiss_data <- read.dta("figure1a.dta")
attach(swiss_data)

pdf(file="C:/Users/ba9ze3-1/Dropbox/Dokumente/Papers/Schweiz IG MP Paper/Final Submission/Giger_Kluever_Replication/figure1a.pdf", height=8, width=11) 
par(mfrow=c(1,1),mar=c(5,5,1.5,3))
hist(econ_no, main = "", col="grey60", xlab = "Average number of sectional groups", ylab = "Number of MPs", 
      cex=1.5, cex.axis=1.5, cex.lab=1.5)
dev.off()


# Figure 1b
# =============================
rm(list=ls(all=TRUE))
setwd("C:/Users/ba9ze3-1/Dropbox/Dokumente/Papers/Schweiz IG MP Paper/Final Submission/Giger_Kluever_Replication")
swiss_data <- read.dta("figure1b.dta")
attach(swiss_data)

pdf(file="C:/Users/ba9ze3-1/Dropbox/Dokumente/Papers/Schweiz IG MP Paper/Final Submission/Giger_Kluever_Replication/figure1b.pdf", height=8, width=11) 
par(mfrow=c(1,1),mar=c(5,5,1.5,3))
hist(ngo_no, main = "", xlim=c(0,20), col="grey60", xlab = "Average number of cause groups", ylab = "Number of MPs", 
      cex=1.5, cex.axis=1.5, cex.lab=1.5)
dev.off()




# =============================
# Figure 2
# =============================

# Figure 2a
# =============================
rm(list=ls(all=TRUE))
setwd("C:/Users/ba9ze3-1/Dropbox/Dokumente/Papers/Schweiz IG MP Paper/Final Submission/Giger_Kluever_Replication")
swiss_data <- read.dta("figure2a.dta")
attach(swiss_data)

pdf(file="C:/Users/ba9ze3-1/Dropbox/Dokumente/Papers/Schweiz IG MP Paper/Final Submission/Giger_Kluever_Replication/figure2a.pdf", height=8, width=11) 
par(mfrow=c(1,1),mar=c(5,5,1.5,3))
barplot(igs_per_party,  ylim=c(0,160), ylab= "Average number of sectional groups", xlab="Political party",  cex=1.5,  cex.axis=1.5, cex.lab=1.5, cex.names=0.7, names.arg=c("FDP", "SVP", "CVP", "SP", "GPS", "FPS", "LPS", "Lega", "LdU",  "-",  "GB", "PdA", "AdG",  "Al", "POCH", "EVP", "Dem", "EDU",  "glp", "CSP", "SD"))  
dev.off()  


# Figure 2b
# =============================
rm(list=ls(all=TRUE))
setwd("C:/Users/ba9ze3-1/Dropbox/Dokumente/Papers/Schweiz IG MP Paper/Final Submission/Giger_Kluever_Replication")
swiss_data <- read.dta("figure2b.dta")
attach(swiss_data)

pdf(file="C:/Users/ba9ze3-1/Dropbox/Dokumente/Papers/Schweiz IG MP Paper/Final Submission/Giger_Kluever_Replication/figure2b.pdf", height=8, width=11) 
par(mfrow=c(1,1),mar=c(5,5,1.5,3))
barplot(igs_per_party, ylim=c(0,160), ylab= "Average number of cause groups", xlab="Political party",  cex=1.5,  cex.axis=1.5, cex.lab=1.5, cex.names=0.7, names.arg=c("SP", "FDP", "CVP", "SVP", "GPS", "EVP", "LPS", "Al", "GB",  "LdU",  "CSP", "glp", "EDU",  "-", "PdA", "FPS", "POCH", "SD",  "Lega", "AdG", "Dem"))  
dev.off()  






# ======================================
# Figure 3
# ======================================

# Figure 3a
# =============================
rm(list=ls(all=TRUE))
setwd("C:/Users/ba9ze3-1/Dropbox/Dokumente/Papers/Schweiz IG MP Paper/Final Submission/Giger_Kluever_Replication")

fixeff_raw <- read.csv("fixeff_m1.csv")
fixeff_untransposed <- fixeff_raw[1,]
fixeff2 <- t(fixeff_untransposed)
fixeff <- fixeff2[1:11]
varcov_raw <- read.csv("varcov_m1.csv")
(varcov <- varcov_raw[1:11,1:11])

# Step 1: get b and V
(b <- as.vector(fixeff))		        
(V <- as.matrix(varcov))		      
(sqrt(diag(V)))                         
                      
# Step 2: Set up MVN distribution
nsim <- 1000
mvn <- mvrnorm(nsim, mu=b, Sigma=V)         
 
# Step 3: Get coefficients
dat <- read.dta("RGraph_Data_m1.dta")
attach(dat)
str(dat)

# Predicted values
pred <- list(rep(NA,31))
v <- seq(0, 30, by=1)
for (i in seq_along(v)){
pred[[i]] <-  (
            + mvn[,1] * i
            + mvn[,2] * mean(dat$ngo_no)
            + mvn[,3] * mean(dat$mp_number)
            + mvn[,4] * median(dat$partei_einig) 
            + mvn[,5] * mean(dat$closevolk)
            + mvn[,6] * mean(dat$closeparl)
            + mvn[,7] * mean(dat$bet)
            + mvn[,8] * mean(dat$months2elec)
            + mvn[,9] * median(dat$RefOblig)
            + mvn[,10] * median(dat$RefFak)
            + mvn[,11] * 1
            )
}

# Transform lists into matrixes
pred_matrix <- do.call(cbind, pred)     
 
# Loop to calculate probabilities
pred_prob <- list(rep(NA,31))
v <- seq(0, 30, by=1)
for (i in seq_along(v)){
  pred_prob[[i]] <- exp(pred_matrix[,i]) / (1 + exp(pred_matrix[,i]))
}

pred_prob_matrix <- do.call(cbind, pred_prob)     

# Taking the means and quantiles of simulated values
pred_prob_mean <- list(rep(NA,31))
pred_prob_lb <- list(rep(NA,31))
pred_prob_ub <- list(rep(NA,31))

v <- seq(0, 30, by=1)
for (i in seq_along(v)){
pred_prob_mean[[i]] <- mean(pred_prob_matrix[,i])
pred_prob_lb[[i]] <- quantile(pred_prob_matrix[,i], c(.05))
pred_prob_ub[[i]] <- quantile(pred_prob_matrix[,i], c(.950))
}


# Plotting predicted values
pdf(file="C:/Users/ba9ze3-1/Dropbox/Dokumente/Papers/Schweiz IG MP Paper/Final Submission/Giger_Kluever_Replication/figure3a.pdf", height=8, width=11)
par(mfrow=c(1,1),mar=c(5,5,1.5,5)) 
hist(dat$econ_no,  xlim=c(0,30), ylab=NULL, xlab=NULL, main=NULL, axes=F, col=rgb(0.95,0.95,0.95, 0.5), lwd=0.1)
mtext("Frequency",side=4,line=4, cex=1.7, las=0) 
axis(4, ylim=c(0,14000), las=0, cex.axis=1.5, cex.lab=1.5)
par(new=TRUE)
plot(1, type="n",  xlim=c(0,30), ylim=c(0,1), xlab= "Number of sectional groups",
      ylab="Predicted probability of defection (incl. 95% CI)", cex=1.7, cex.axis=1.5, cex.lab=1.5)
      lines(0:30,pred_prob_mean, col="black", lwd=2)
      lines(0:30,pred_prob_ub, col="black", lty=2, lwd=4)
      lines(0:30,pred_prob_lb, col="black", lty=2, lwd=4)
dev.off()





# Figure 3b
# =============================

rm(list=ls(all=TRUE))
setwd("C:/Users/ba9ze3-1/Dropbox/Dokumente/Papers/Schweiz IG MP Paper/Final Submission/Giger_Kluever_Replication/")

fixeff_raw <- read.csv("fixeff_m1.csv")
fixeff_untransposed <- fixeff_raw[1,]
fixeff2 <- t(fixeff_untransposed)
fixeff <- fixeff2[1:11]
varcov_raw <- read.csv("varcov_m1.csv")
(varcov <- varcov_raw[1:11,1:11])

# Step 1: get b and V
(b <- as.vector(fixeff))		        
(V <- as.matrix(varcov))		      
(sqrt(diag(V)))                         
                      
# Step 2: Set up MVN distribution
nsim <- 1000
mvn <- mvrnorm(nsim, mu=b, Sigma=V)         

# Step 3: Get coefficients
dat <- read.dta("RGraph_Data_m1.dta")
attach(dat)
str(dat)

# Predicted values
pred <- list(rep(NA,31))
v <- seq(0, 30, by=1)
for (i in seq_along(v)){
pred[[i]] <-  (
            + mvn[,1] * mean(dat$econ_no)
            + mvn[,2] * i
            + mvn[,3] * mean(dat$mp_number)
            + mvn[,4] * median(dat$partei_einig) 
            + mvn[,5] * mean(dat$closevolk)
            + mvn[,6] * mean(dat$closeparl)
            + mvn[,7] * mean(dat$bet)
            + mvn[,8] * mean(dat$months2elec)
            + mvn[,9] * median(dat$RefOblig)
            + mvn[,10] * median(dat$RefFak)
            + mvn[,11] * 1
            )
}

# Transform lists into matrixes
pred_matrix <- do.call(cbind, pred)     

# Loop to calculate probabilities
pred_prob <- list(rep(NA,31))
v <- seq(0, 30, by=1)
for (i in seq_along(v)){
  pred_prob[[i]] <- exp(pred_matrix[,i]) / (1 + exp(pred_matrix[,i]))
}

pred_prob_matrix <- do.call(cbind, pred_prob)     

# Taking the means and quantiles of simulated values
pred_prob_mean <- list(rep(NA,31))
pred_prob_lb <- list(rep(NA,31))
pred_prob_ub <- list(rep(NA,31))

v <- seq(0, 30, by=1)
for (i in seq_along(v)){
pred_prob_mean[[i]] <- mean(pred_prob_matrix[,i])
pred_prob_lb[[i]] <- quantile(pred_prob_matrix[,i], c(.025))
pred_prob_ub[[i]] <- quantile(pred_prob_matrix[,i], c(.975))
}

# Plotting predicted values
pdf(file="C:/Users/ba9ze3-1/Dropbox/Dokumente/Papers/Schweiz IG MP Paper/Final Submission/Giger_Kluever_Replication/figure3b.pdf", height=8, width=11)
par(mfrow=c(1,1),mar=c(5,5,1.5,5)) 
hist(dat$ngo_no,  xlim=c(0,30), ylab=NULL, xlab=NULL, main=NULL, axes=F, col=rgb(0.95,0.95,0.95, 0.5), lwd=0.1)
mtext("Frequency",side=4,line=4, cex=1.7, las=0) 
axis(4, ylim=c(0,14000), las=0, cex.axis=1.5, cex.lab=1.5)
par(new=TRUE)
plot(1, type="n",  xlim=c(0,30), ylim=c(0,1), xlab= "Number of cause groups",
      ylab="Predicted probability of defection (incl. 95% CI)", cex=1.7, cex.axis=1.5, cex.lab=1.5)
      lines(0:30,pred_prob_mean, col="black", lwd=2)
      lines(0:30,pred_prob_ub, col="black", lty=2, lwd=4)
      lines(0:30,pred_prob_lb, col="black", lty=2, lwd=4)
dev.off()







# ======================================
# Figure 4
# ======================================

# Figure 4a
# ==============================
rm(list=ls(all=TRUE))
setwd("C:/Users/ba9ze3-1/Dropbox/Dokumente/Papers/Schweiz IG MP Paper/Final Submission/Giger_Kluever_Replication")

fixeff_raw <- read.csv("fixeff_m2.csv")
fixeff_untransposed <- fixeff_raw[1,]
fixeff2 <- t(fixeff_untransposed)
fixeff <- fixeff2[1:11]
varcov_raw <- read.csv("varcov_m2.csv")
(varcov <- varcov_raw[1:11,1:11])

# Step 1: get b and V
(b <- as.vector(fixeff))		        
(V <- as.matrix(varcov))		      
(sqrt(diag(V)))                         
                      
# Step 2: Set up MVN distribution
nsim <- 1000
mvn <- mvrnorm(nsim, mu=b, Sigma=V)         

# Step 3: Get coefficients
dat <- read.dta("RGraph_Data_m2.dta")
attach(dat)
str(dat)

# Predicted values
pred <- list(rep(NA,251))
v <- seq(0, 250, by=1)
for (i in seq_along(v)){
pred[[i]] <-  (
            + mvn[,1] * i
            + mvn[,2] * mean(dat$party_ngo_no)
            + mvn[,3] * mean(dat$mp_number)
            + mvn[,4] * median(dat$partei_einig) 
            + mvn[,5] * mean(dat$closevolk)
            + mvn[,6] * mean(dat$closeparl)
            + mvn[,7] * mean(dat$bet)
            + mvn[,8] * mean(dat$months2elec)
            + mvn[,9] * median(dat$RefOblig)
            + mvn[,10] * median(dat$RefFak)
            + mvn[,11] * 1
            )
}

# Transform lists into matrixes
pred_matrix <- do.call(cbind, pred)     

# Loop to calculate probabilities
pred_prob <- list(rep(NA,251))
v <- seq(0, 250, by=1)
for (i in seq_along(v)){
  pred_prob[[i]] <- exp(pred_matrix[,i]) / (1 + exp(pred_matrix[,i]))
}

pred_prob_matrix <- do.call(cbind, pred_prob)     

# Taking the means and quantiles of simulated values
pred_prob_mean <- list(rep(NA,251))
pred_prob_lb <- list(rep(NA,251))
pred_prob_ub <- list(rep(NA,251))

v <- seq(0, 250, by=1)
for (i in seq_along(v)){
pred_prob_mean[[i]] <- mean(pred_prob_matrix[,i])
pred_prob_lb[[i]] <- quantile(pred_prob_matrix[,i], c(.025))
pred_prob_ub[[i]] <- quantile(pred_prob_matrix[,i], c(.975))
}

# Plotting predicted values
pdf(file="C:/Users/ba9ze3-1/Dropbox/Dokumente/Papers/Schweiz IG MP Paper/Final Submission/Giger_Kluever_Replication/figure4a.pdf", height=8, width=11)
par(mfrow=c(1,1),mar=c(5,5,1.5,5)) 
hist(dat$party_econ_no,  xlim=c(0,250), ylab=NULL, xlab=NULL, main=NULL, axes=F, col=rgb(0.95,0.95,0.95, 0.5), lwd=0.1)
mtext("Frequency",side=4,line=4, cex=1.7, las=0) 
axis(4, ylim=c(0,4000), las=0, cex.axis=1.5, cex.lab=1.5)
par(new=TRUE)
plot(1, type="n",  xlim=c(0,250), ylim=c(0,1), xlab= "Number of sectional groups",
      ylab="Predicted probability of defection (incl. 95% CI)", cex=1.7, cex.axis=1.5, cex.lab=1.5)
      lines(0:250,pred_prob_mean, col="black", lwd=2)
      lines(0:250,pred_prob_ub, col="black", lty=2, lwd=4)
      lines(0:250,pred_prob_lb, col="black", lty=2, lwd=4)
dev.off()






# Figure 4b
# ==============================
rm(list=ls(all=TRUE))
setwd("C:/Users/ba9ze3-1/Dropbox/Dokumente/Papers/Schweiz IG MP Paper/Final Submission/Giger_Kluever_Replication/")

fixeff_raw <- read.csv("fixeff_m2.csv")
fixeff_untransposed <- fixeff_raw[1,]
fixeff2 <- t(fixeff_untransposed)
fixeff <- fixeff2[1:11]
varcov_raw <- read.csv("varcov_m2.csv")
(varcov <- varcov_raw[1:11,1:11])

# Step 1: get b and V
(b <- as.vector(fixeff))		        
(V <- as.matrix(varcov))		      
(sqrt(diag(V)))                         
                      
# Step 2: Set up MVN distribution
nsim <- 1000
mvn <- mvrnorm(nsim, mu=b, Sigma=V)         

# Step 3: Get coefficients
dat <- read.dta("RGraph_Data_m2.dta")
attach(dat)
str(dat)

# Predicted values
pred <- list(rep(NA,251))
v <- seq(0, 250, by=1)
for (i in seq_along(v)){
pred[[i]] <-  (
            + mvn[,1] * mean(dat$party_econ_no)
            + mvn[,2] * i
            + mvn[,3] * mean(dat$mp_number)
            + mvn[,4] * median(dat$partei_einig) 
            + mvn[,5] * mean(dat$closevolk)
            + mvn[,6] * mean(dat$closeparl)
            + mvn[,7] * mean(dat$bet)
            + mvn[,8] * mean(dat$months2elec)
            + mvn[,9] * median(dat$RefOblig)
            + mvn[,10] * median(dat$RefFak)
            + mvn[,11] * 1
            )
}

# Transform lists into matrixes
pred_matrix <- do.call(cbind, pred)     

# Loop to calculate probabilities
pred_prob <- list(rep(NA,251))
v <- seq(0, 250, by=1)
for (i in seq_along(v)){
  pred_prob[[i]] <- exp(pred_matrix[,i]) / (1 + exp(pred_matrix[,i]))
}

pred_prob_matrix <- do.call(cbind, pred_prob)     

# Taking the means and quantiles of simulated values
pred_prob_mean <- list(rep(NA,251))
pred_prob_lb <- list(rep(NA,251))
pred_prob_ub <- list(rep(NA,251))

v <- seq(0, 250, by=1)
for (i in seq_along(v)){
pred_prob_mean[[i]] <- mean(pred_prob_matrix[,i])
pred_prob_lb[[i]] <- quantile(pred_prob_matrix[,i], c(.025))
pred_prob_ub[[i]] <- quantile(pred_prob_matrix[,i], c(.975))
}


# Plotting predicted values
pdf(file="C:/Users/ba9ze3-1/Dropbox/Dokumente/Papers/Schweiz IG MP Paper/Final Submission/Giger_Kluever_Replication/figure4b.pdf", height=8, width=11)
par(mfrow=c(1,1),mar=c(5,5,1.5,5)) 
hist(dat$party_ngo_no,  xlim=c(0,250), ylab=NULL, xlab=NULL, main=NULL, axes=F, col=rgb(0.95,0.95,0.95, 0.5), lwd=0.1)
mtext("Frequency",side=4,line=4, cex=1.7, las=0) 
axis(4, ylim=c(0,4000), las=0, cex.axis=1.5, cex.lab=1.5)
par(new=TRUE)
plot(1, type="n",  xlim=c(0,250), ylim=c(0,1), xlab= "Number of cause groups",
      ylab="Predicted probability of defection (incl. 95% CI)", cex=1.7, cex.axis=1.5, cex.lab=1.5)
      lines(0:250,pred_prob_mean, col="black", lwd=4)
      lines(0:250,pred_prob_ub, col="black", lty=2, lwd=4)
      lines(0:250,pred_prob_lb, col="black", lty=2, lwd=4)
dev.off()

